

clear all; close all; dbclear if error; 

%% ESTIMATION OPTIONS
% quarterly estimation, all cities, locations only
options.SaveFigures  = 'No' ;   % Yes or No
options.SaveEstimates = 'No' ; % Yes or No
 
options.RunDelta = 0 ;
options.BetaMax = 0.3 ;
options.Na = 50 ;
options.learning = 'Learning' ; % 'Learning' or 'NoLearning' ;

options.HousePriceElasticity = 1 ;

%% COLOR SCHEME
MyOrange = [ 231,101,45 ] ;
MyOrange = MyOrange / 255 ;
MyBlue =  [ 60,106,209 ] ;
MyBlue = MyBlue / 255 ;


%% LOAD DATA

% migration rate
data1   = csvread('source/QuarterlyMigProba.csv'); 
data1_1 = csvread('source/LFExitRate.csv') ;
data2   = csvread('source/LaborMarketBasics.csv');
data3   = csvread('source/LaborMarketBasics2periods.csv');
data4   = csvread('source/IndEmpIV.csv'); 
data5   = csvread('source/MigShares.csv'); 
data6   = csvread('source/HousePricesArea.csv'); 
data10  = csvread('source/EUyear1.csv'); 
data11  = csvread('source/ANPE.csv');
data12  = csvread('source/WageTenure.csv'); 

% Migration rate
MigQ = data1(1) ;

% labor force exit rate
LFexitQ = data1_1(1) ;

% EU and UE data: QUARTERLY FREQUENCY
City    = data2(:,1) ;
EU      = data2(:,2) ;
UE      = data2(:,3) ;
W       = data2(:,4) ;
wt      = data2(:,5) ;
u       = data2(:,6) ;

% split sample in two subperiods
u0   = data3(:,2);
U0   = data3(:,3);
EU0  = data3(:,4);
UE0  = data3(:,5);
W0   = data3(:,6);
pop0 = data3(:,7); pop0 = pop0 ./ sum( pop0 ) ;
u1   = data3(:,8);
U1   = data3(:,9);
EU1  = data3(:,10);
UE1  = data3(:,11);
W1   = data3(:,12);
pop1 = data3(:,13); pop1 = pop1 ./ sum( pop1 ) ;

% Bartik industry employment IV
indemp0 = data4(:,2) ;
indemp1 = data4(:,3) ;

% migration shares
mig0    = data5(:,2) ;
mig1    = data5(:,4) ;
 
% House prices and amenities for over-id
area       = data6(:,2) ;
houseprice = data6(:,3) ; % from notaires

% dhat data: ANNUAL FREQUENCY
EU1Y    = data10(:,3) ;
EUmY    = data10(:,4) ;

% ANPE and search data
FractionANPE            = data11(1);
UnempDuration           = data11(2) ;
QuartersSinceContact    = data11(3) ;
QuartersSinceOffer      = data11(4) ;

% Wages by tenure data
domegaData     = data12(:,2) ; 
dSumOmegaData  = data12(:,3) ; 

% Popolation weights : density
pop = wt .* ( 1 + EU ./ ( EU + UE ) ) ;
pop = pop / sum(pop) ; 


%% PRE-SET PARAMETERS FROM OFFLINE ESTIMATION AND NORMALIZATION
Params.phi = 0.1 ; 
Params.psi = Params.phi / ( 1 - Params.phi ) ;
Params.omega = 0.23 ;
Params.b = 1 ;
Params.phiLearning = 0.01 ;

% Correct migration proba for time aggregation
Params.mu = - log( 1 - MigQ ) ;
muA       = 4 * Params.mu ;

% Effective discount rate in adjusted Bellman equation
rhoAnnual = 0.05 ;
Params.rho = rhoAnnual / 4 + Params.mu + Params.phiLearning ;
   

%% ADJUST FOR TIME AGGREGATION

% Quarterly instantaneous rates

FpSQ = - log( 1 - EU - UE ) ;
SQ = FpSQ ./ ( EU + UE ) .* EU ;
FQ = FpSQ ./ ( EU + UE ) .* UE ;

FpSQ0 = - log( 1 - EU0 - UE0 ) ;
SQ0 = FpSQ0 ./ ( EU0 + UE0 ) .* EU0 ;
FQ0 = FpSQ0 ./ ( EU0 + UE0 ) .* UE0 ;

FpSQ1 = - log( 1 - EU1 - UE1 ) ;
SQ1 = FpSQ1 ./ ( EU1 + UE1 ) .* EU1 ;
FQ1 = FpSQ1 ./ ( EU1 + UE1 ) .* UE1 ;

s = SQ ;
f = FQ ;

s0 = SQ0 ;
s1 = SQ1 ;
f0 = FQ0 ;
f1 = FQ1 ;
    
st = s - Params.mu ;
st0 = s0 - Params.mu ;
st1 = s1 - Params.mu ;

% annual instantaneous rates
SA = 4 * SQ ;
FA = 4 * FQ ;
SAt = SA - muA ;


%% ESTIMATION OF LABOR FORCE EXIT RATE (Delta)

% Find Delta from labor force exit rate
Params.Delta = LFexitQ ; 

%% ESTIMATION OF DELTA HAT

% Adjust for zeros
EU1Y(EU1Y == 0 ) = min( EU1Y(EU1Y > 0 ) ) ;

% Grid for delta hat
Ng = 200 ;
dhatMin = 0.2 ; 
dhatMax = 3 ;
dhatGrid = linspace(dhatMin,dhatMax,Ng)'; 
predictedEU1Y  = zeros(size(SQ,1),Ng) ;
EU1Mat  = repmat( EU1Y , 1 , Ng ) - muA ;

for i=1:Ng
    
    d = dhatGrid(i) ;
    D = Params.Delta + Params.mu ;
    S = EUmY - muA ;
    Sd = S / d ;
    sq = sqrt( 2*D + d^2 ) ;
    s2 = sqrt(2) ;
    
    prefactor = Sd ./ ( 2 * D + 2 * S - Sd.^2 ) ;
    term1     = sq * erf( sq / s2 ) ;
    term2     = ( d - Sd ) ...
             .* ( 1 - exp( Sd.^2/2 - D - S ) ...
                   .* ( 1 + erf( ( d - Sd ) / s2 ) ) ...
                ) ;
   
    pred = prefactor .* ( term1 + term2 ) ;
     
    predictedEU1Y(:,i) = pred ;
    
end

objMat = ( predictedEU1Y - EU1Mat ) ./ EU1Mat ;
objGrid = sum( objMat.^2 .* repmat(pop,1,Ng) , 1 ) ;

% Old stuff
SAMat   = repmat( EUmY - muA , 1 , Ng ) ;
dhatMat = repmat( dhatGrid' , size(SAt,1) , 1 ) ;
xMat = SAMat ./ dhatMat ;
s2 = sqrt(2) ;

survivors = 1 ./ ( 2 * dhatMat - xMat ) ...
         .* (   dhatMat .* erfc( dhatMat / s2 ) ...
              + ( xMat - dhatMat) .* exp( xMat.^2/2 - SAMat ) ...
                .* ( erfc( ( dhatMat - xMat ) / s2 ) - 2  )  ...
            ) ;


objMatOld = ( 1 - survivors - EU1Mat ) ./ EU1Mat ;
objGridOld = sum( objMatOld.^2 .* repmat(pop,1,Ng) , 1 ) ;
[~,iminOld] = min(objGridOld);
dhatAnnualOld = dhatGrid(iminOld);

figure()
set(gcf,'Position',[100 100 500 500])
fontsize = 14 ;
plot( dhatGrid / 2 , objGrid  , 'LineWidth',2); hold on;
xlim([0.2,1])
xlabel('$\hat{\delta}$','interpreter','latex')
ylabel('sum of square residuals')
title('Identification of $\hat{\delta}$','interpreter','latex')
AxisFonts(0.75*fontsize,fontsize,0.75*fontsize,fontsize)
hold off

[~,imin] = min(objGrid);
dhatAnnual = dhatGrid(imin);

% Predicted
EU1pred = predictedEU1Y(:,imin) ;

figure()
set(gcf,'Position',[100 100 500 480])
    
fontsize = 20 ;
maxplot = max(max(EU1Y,EU1pred)) ;
minplot = min(min(EU1Y,EU1pred)) ;
scatter(EU1Y - muA , EU1pred , 5 * 1e4 * pop ,  MyBlue , 'LineWidth',1.5); hold on;
plot([ 0.8 * minplot ; 1.05 * maxplot ] ,[ 0.8 * minplot ; 1.05 * maxplot ],'Color',MyOrange, 'LineWidth',2)
xlim([ 0.8 * minplot ; 0.8 * maxplot ])
ylim([ 0.8 * minplot ; 0.8 * maxplot ])
xlabel('Data')
ylabel('Model')
grid on
xlim([0.04,0.15])
ylim([0.04,0.15])
xh = get(gca,'xlabel') ;
px = get(xh,'position') ;
px(2) = px(2) - 0.002 ;
set(xh,'position',px) ;
yh = get(gca,'ylabel') ;
py = get(yh,'position') ;
py(1) = py(1) - 0.002 ;
set(yh,'position',py) ;
AxisFonts(0.95*fontsize,fontsize,0.95*fontsize,fontsize)

if strcmp(options.SaveFigures,'Yes')==1
    
    % pdf
    filename = 'results/JobLossRatesFirstYear' ;
    run printPDF.m ;
    
    % tiff
    exportgraphics(gcf,'results/JobLossRatesFirstYear.tiff','Resolution',1000) 
    
end

hold off


% Adjust dhat for annual vs. quarterly frequency
dhat = dhatAnnual / 2 ;
Params.dhat = dhat ;

fprintf('\n')
disp('********')
fprintf('\n')
disp(['Estimate of dhat is = ' num2str(dhat,4)])
fprintf('\n')
disp('********')


%% ESTIMATE DELTA

Nt   = 25 ; % number of points: years w/ 0
tMax = 4 * 24 ; % number of points: years
tGrid = (0:4:tMax)' ; tGrid(1) = 0.5 ;
Ns = size(s,1) ;

Mass       = zeros(Ns,Nt) ;
SumW       = zeros(Ns,Nt) ;
MeanW      = zeros(Ns,Nt) ;

D = Params.mu + Params.Delta ;

for it=1:Nt

    t = tGrid(it) ;

    for i=1:Ns

        disp(['(t,i) = (' num2str(t) ',' num2str(i) ')'])

        ss = st(i) ;

        g = @(z) dhat * exp( - dhat^2 / ( 2 * t ) * ( z + t ).^2 ) / sqrt( 2 * pi * t ) ;

        h0G = @(y,y0) ...
              ss * exp( - ss * y0 ) ...
           .* exp( - D * t ) .* ( g( y - y0) - exp( 2 * dhat^2 * y0 ) .* g( y + y0 ) ) ;

        h0Gyy0 = @(y,y0) h0G(y,y0) .* (y - y0) ;

        % compute integrals
        bound = 100 ;
        
        Mass(i,it)  = integral2( h0G    , 0 , 2 * bound , 0 , bound ) ;
        SumW(i,it)  = integral2( h0Gyy0 , 0 , 2 * bound , 0 , bound ) ;
        
        % compute means
        MeanW(i,it) = SumW(i,it) / Mass(i,it) ;    

    end


end
    
omega = sum( MeanW .* repmat( pop , [ 1 Nt ] ) , 1 )' ;
domega = omega - omega(1) ;

Nd = 1000 ;
c1Grid = linspace(0,0.02,Nd)' ;
obj = zeros(Nd,1) ;

for i=1:Nd
    
    obj(i) = sum( ( domegaData - c1Grid(i) * domega ).^2 ) / size(domega,1) ;
    
end



figure()
set(gcf,'Position',[100 100 500 500])
fontsize = 14 ;
plot(c1Grid,obj,'LineWidth',2)
xlim([0,0.02])
xlabel('$\frac{\delta}{(1-\beta)\rho + \beta}$','interpreter','latex')
ylabel('sum of square residuals')
AxisFonts(0.75*fontsize,fontsize,0.75*fontsize,fontsize)
hold off

[~,ic1] = min(obj) ;
c1 = c1Grid(ic1) ;



figure()
set(gcf,'Position',[100 100 500 500])
fontsize = 14 ;
plot(tGrid'/4,domegaData,'Color',MyBlue); hold on;
plot(tGrid'/4,c1 * domega,'--','Color',MyOrange); hold on;
legend('Data','Model','Fontsize',fontsize)
legend boxoff
ylim([-0.8 0.2])

xlabel('Years at job')
ylabel('Within-spell wages relative to new job wages')
title('Identification of $\delta$','interpreter','latex')
AxisFonts(0.75*fontsize,fontsize,0.75*fontsize,fontsize)
hold off


%% ESTIMATE BETA

Targets.ls = 0.37 ; % using wages net of "cotisations patronales" in FICUS data

Nb = 1000 ;
betaGrid = linspace(0.01,0.99,Nb)' ;
obj = zeros(Nb,1) ;
LaborShare = zeros(Nb,1) ;

for i=1:Nb
    
    beta = betaGrid(i) ;
    delta = c1 * ( ( 1 - beta ) * Params.rho + beta ) / beta ;
    sigma = delta / dhat ;
    
    kappa = ( dhat + sqrt( D + dhat^2) ) / sigma ;
    tau = delta / sigma^2 * ( sqrt(1 + 2 * Params.rho * sigma^2 / delta^2 ) - 1 ) ; 
    zl0   = tau * ( Params.rho + delta - sigma^2/2) / (1 + tau) ;

    zeta = st / delta ;
    Ez = zl0 ./ Params.rho .* kappa .* zeta ./ ( zeta - 1 ) ./ ( kappa - 1 ) ;

    beta = betaGrid(i) ;
    LaborShare(i) = beta + ( 1 - beta) * sum( pop ./ Ez ) ;
    obj(i) = abs( LaborShare(i) - Targets.ls ) / Targets.ls ;
    
end


figure()
set(gcf,'Position',[100 100 500 500])
fontsize = 14 ;
plot(betaGrid,LaborShare - Targets.ls,'LineWidth',2)
xlim([0.025,0.2])
xlabel('$\beta$','interpreter','latex')
ylabel('model-data labor share difference')
title('Identification of $\beta$','interpreter','latex')
AxisFonts(0.75*fontsize,fontsize,0.75*fontsize,fontsize)
hold off

[~,ib] = min(obj) ;
Params.beta = betaGrid(ib) ;
Params.delta = c1 * ( ( 1 - Params.beta ) * Params.rho + Params.beta ) ;
Params.sigma = Params.delta / dhat ;

fprintf('\n')
disp('********')
fprintf('\n')
disp(['Estimate of beta is = ' num2str(Params.beta,4)])
disp(['Estimate of delta is = ' num2str(Params.delta,4)])
disp(['Estimate of sigma is = ' num2str(Params.sigma,4)])
fprintf('\n')
disp('********')

%% DEFINE REMAINING TARGETS AND USEFUL DATA

Targets.UnempDuration = UnempDuration ;
Targets.FractionANPE = FractionANPE ; 
Targets.QuartersSinceContact = QuartersSinceContact ;
Targets.QuartersSinceOffer = QuartersSinceOffer ;

data.s = s ;
data.st = st ;
data.f = f ;
data.u = u ;
data.W = W ;
data.pop = pop ;

data.mig0 = mig0 ;
data.mig1 = mig1 ;
data.indemp0 = indemp0 ;
data.indemp1 = indemp1 ;

data.u0 = u0 ;
data.s0 = s0 ;
data.st0 = st0 ;
data.f0 = f0 ;
data.W0 = W0 ;
data.pop0 = pop0 ;

data.u1 = u1 ;
data.s1 = s1 ;
data.st1 = st1 ;
data.f1 = f1 ;
data.W1 = W1 ;
data.pop1 = pop1 ;

data.houseprice = houseprice ; 
data.area = area ;


%% ESTIMATE OTHER PARAMETERS : NO HOSIOS CONDITION (DECENTRALIZED EQUILIBRIUM)

options.EnforceHosios = 0 ;

[ ParamsNH , InversionNH ] = EstimateRest( options , Params , Targets , data ) ;
ParamsNH

%% ESTIMATE OTHER PARAMETERS : HOSIOS CONDITION + NO EXTERNALITY (PLANNER ALLOCATION/DIRECTED SEARCH)

options.EnforceHosios = 1 ;

Na = options.Na ;
alphaGrid = linspace( 0.01 , 0.7 , Na )' ;

alphaUpdateUW = zeros( Na , 1 ) ;
gammaUpdateUW = zeros( Na , 1 ) ;

% Loop over regressions
for i = 1:Na
    
    disp([ 'i = ' num2str(i) ] )
    
    Params.beta = alphaGrid(i) ;
    options.betamax = 0.6 ;
    
    [ ParamsHosios , InversionHosios ] = EstimateRest( options , Params , Targets , data ) ;
    ParamsStore{i} = ParamsHosios ;
    InversionStore{i} = InversionHosios ;
    
    SbarHosios = @(zzeta) ParamsHosios.tau ./ ( zzeta - 1 ) ./ ( ParamsHosios.tau + zzeta ) ;

    % Prepare input to recover relevant variables for alpha regression
    Data0.u     = u0 ;
    Data0.s     = s0 ;
    Data0.st    = st0 ;
    Data0.f     = f0 ;
    Data0.W     = W0 ;
    Data0.pop   = pop0 ;
    Data0.area  = area ;

    Data1.u     = u1 ;
    Data1.s     = s1 ;
    Data1.st    = st1 ;
    Data1.f     = f1 ;
    Data1.W     = W1 ;
    Data1.pop   = pop1 ;
    Data1.area  = area ;
    
    model0 = inversion( Data0 , ParamsHosios , options ) ;
    model1 = inversion( Data1 , ParamsHosios , options ) ;
    
    Data0.HousePrice = model0.HousePrice ;
    Data0.lL = model0.lL ;
    Data0.Ez = model0.Ez ;
    
    Data1.HousePrice = model1.HousePrice ;
    Data1.lL = model1.lL ;
    Data1.Ez = model1.Ez ;
    
    [ ~ , lX0 ] = inversionAm( Data0 , ParamsHosios  ) ;
    [ ~ , lX1 ] = inversionAm( Data1 , ParamsHosios  ) ;
    
    P2 = ( ParamsHosios.omega + ParamsHosios.eps ) / ParamsHosios.P1 ;
    
    % Prepare regression variables
    
    fzeta0 = (model0.zetaR).^( ParamsHosios.p1zeta - 1 ) .* ( 1 - model0.zetaR ).^( ParamsHosios.p2zeta - 1 ) ;
    fzeta1 = (model1.zetaR).^( ParamsHosios.p1zeta - 1 ) .* ( 1 - model1.zetaR ).^( ParamsHosios.p2zeta - 1 ) ;
    fx0    = exp( - lX0.^2 / 2 / ParamsHosios.stdx^2 ) ;
    fx1    = exp( - lX1.^2 / 2 / ParamsHosios.stdx^2 ) ;
    
    dzeta0 = ones( size(s0) ) ;
    dzeta1 = ones( size(s1) ) ;
    
    LHS0  = log( f0 .* (ParamsHosios.B0).^( - model0.zeta ) .* ( (model0.bpv).^(model0.zeta) ) ) ;
    LHS1  = log( f1 .* (ParamsHosios.B0).^( - model0.zeta ) .* ( (model1.bpv).^(model1.zeta) ) ) ;
    
    RHSa0  = log( fzeta0 .* dzeta0 ./ model0.Um ./ fx0 ) ;
    RHSa1  = log( fzeta1 .* dzeta1 ./ model1.Um ./ fx1 ) ;
    
    RHSg0 = log(    ( ParamsHosios.B0).^(model0.zeta) ...
                 .* (model0.bpv).^( 1 - model0.zeta - ParamsHosios.psi / ParamsHosios.P1 )  ...
     			 .*	(model0.G).^( - ParamsHosios.eps * ParamsHosios.psi / ParamsHosios.P1 )  ...
                 .* SbarHosios( model0.zeta ) .* exp( lX0 + P2 * model0.lH ) ) ;
    RHSg1 = log(    ( ParamsHosios.B0).^(model1.zeta) ...
                 .* (model1.bpv).^( 1 - model1.zeta - ParamsHosios.psi / ParamsHosios.P1 )  ...
     			 .*	(model1.G).^( - ParamsHosios.eps * ParamsHosios.psi / ParamsHosios.P1 )  ...
                 .* SbarHosios( model1.zeta ) .* exp( lX1 + P2 * model1.lH  ) ) ;

    IV0 = log( indemp0 ) ;
    IV1 = log( indemp1 ) ;
    
    % Define first difference
    dLHS  = LHS1 - LHS0 ;
    dRHSa = RHSa1 - RHSa0 ;
    dRHSg = RHSg1 - RHSg0 ;
    dIV   = IV1 - IV0 ;
    
    % Clean from extreme values
    drop = (   isnan( dLHS ) + isnan( dRHSa ) + isnan( dRHSg ) + isnan( dIV ) ...
             + isinf( dLHS ) + isinf( dRHSa ) + isinf( dRHSg ) + isinf( dIV ) ) > 0 ;

    dLHS(drop) = [] ;
    dRHSa(drop) = [] ;
    dRHSg(drop) = [] ;
    dIV(drop) = [] ;
    
    % More IVs
    dIV2 = dIV > mean( dIV ) ;
    
    % REGRESSIONS
    TSLSUW = MyTsls( dLHS , [ dRHSa , dRHSg ] , [ dIV , dIV2 ] , ones(size(dLHS,1)) ) ;
    [ alpha , gamma ] = RecoverAlphaGamma( TSLSUW.Coefficients.Estimate(2) , ...
                                           TSLSUW.Coefficients.Estimate(3) ) ;
    alphaUpdateUW(i) = alpha ;
    gammaUpdateUW(i) = gamma ;

    
end


figure()
nx = 1 ; ny = 2 ;

subplot(nx,ny,1)
plot(alphaGrid,alphaUpdateUW,'LineWidth',2); hold on;
plot(alphaGrid,alphaGrid,'k','LineWidth',1)
ylim([ 0.01 , 1 ] )
xlabel('Initial guess for bargaining power','interpreter','latex')
ylabel('Implied $\alpha$','interpreter','latex')
title('Identification of $\alpha$','interpreter','latex')
hold off

subplot(nx,ny,2)
plot(alphaGrid,gammaUpdateUW,'LineWidth',2); hold on;
xlabel('Initial guess for bargaining power','interpreter','latex')
ylabel('Implied $\gamma$','interpreter','latex')
title('Identification of $\gamma$','interpreter','latex')
ylim( [ 0 , 15 ] )
xlim([0,0.4] )
hold off

hold off

[~ , iL , ~ , ~ ] = FindZeroInterp( alphaGrid - alphaUpdateUW , alphaGrid ) ;

ParamsH = ParamsStore{iL} ; 
InversionH = InversionStore{iL} ; 
ParamsH.alpha = alphaUpdateUW(iL) ;
ParamsH.gamma = gammaUpdateUW(iL) ;

% Get alpha for the non-hosios case by taking the alpha estimate for alphaGrid = beta
[~ , iLNH , ~ , ~ ] = FindZeroInterp( ParamsNH.beta - alphaGrid , alphaGrid ) ;
ParamsNH.alpha = alphaUpdateUW(iLNH) ;
ParamsNH.gamma = gammaUpdateUW(iLNH) ;

disp(['DE: alpha = ' num2str(ParamsNH.alpha) ])
disp(['DE: gamma = ' num2str(ParamsNH.gamma) ])
disp(['SP: alpha = ' num2str(ParamsH.alpha) ])
disp(['SP: gamma = ' num2str(ParamsH.gamma) ])

if strcmp(options.SaveEstimates,"Yes")

    save('estimates/ParamsNH.mat','ParamsNH')
    save('estimates/ParamsH.mat','ParamsH')

    disp(' ')
    disp('*** Saved estimates ***')
    disp(' ')

end


%% OVER-ID FIGURE FOR HOUSING PRICES

% Get data
weights = data.pop ; 

% House prices
keepit = data.houseprice > 0 ;
modelLogPrice = log( InversionNH.HousePrice(keepit) ) ; 
dataLogPrice = log( data.houseprice(keepit) ) ;
wt = weights(keepit) ;
wt = wt / sum(wt) ;

% Normalize mean
modelLogPrice = modelLogPrice - sum( wt .* modelLogPrice ) ;
dataLogPrice  = dataLogPrice  - sum( wt .* dataLogPrice ) ;


figure()
set(gcf,'Position',[100 100 500 480])
fontsize = 20 ;
scatter( dataLogPrice , modelLogPrice , ...
         2*1e4 * wt , ...
         MyBlue ,  'LineWidth', 1.5 ); hold on
plot([ -1 ; 1 ],[ -1 ; 1 ],'LineStyle','-','Color',MyOrange,'LineWidth',2) 
xlabel('Data')
ylabel('Model')
xlim([ -1 ; 1 ])
ylim([ -1 ; 1 ])
AxisFonts(0.95*fontsize,fontsize,0.95*fontsize,fontsize)
grid on

set(gca,'box','off')

if strcmp(options.SaveFigures,'Yes')==1
    
    % pdf
    filename = 'results/FigureHousePrices' ;
    run printPDF.m ;
    
    % tiff
    exportgraphics(gcf,'results/FigureHousePrices.tiff','Resolution',1000) 
    
end

hold off


%% WRITE RESULTS TO DISK

% write table of estimates
load('estimates/ParamsNH.mat'); 
WriteEstimatesTable( ParamsNH , 'Estimates' )

% Export amenities
csvwrite('results/Amenities.csv', [ City , InversionNH.lA ] )









